Two dimensional modulational instability in photorefractive media 
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We study theoretically and experimentally the modulational instability of broad optical beams in 
photorefractive nonlinear media. We demonstrate the impact of the anisotropy of the nonlinearity on 
the growth rate of periodic perturbations. Our findings are confirmed by experimental measurements 
in a strontium barium niobate photorefractive crystal. 
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A plane wave propagating in a medium with focusing 
nonlinearity is unstable with respect to the generation 
of small scale filaments ilj. This so called modulational 
instability (MI) phenomenon, has been extensively stud- 
ied because of its importance as a factor limiting the 
propagation of high power beams. Filament ation may 
also be identified as the first stage in the development of 
turbulent fluctuations in the transverse profile of a laser 
beam 2]. In addition, MI is often considered as a pre- 
cursor for the formation of spatial and/or temporal op- 
tical solitons. As far as optics is concerned, MI has been 
studied in media with various mechanisms of nonlinear 
response including cubic Q , quadratic 3] , nonlocal 0, 
and inertial 6, 7] types of nonlinearity. Importantly, MI 
is not restricted to nonlinear optics but has also been 
studied in many other nonlinear systems including fluids 
[S], plasmas 9] and matter waves 10]. 

In the context of optical beam propagation in nonlin- 
ear materials MI has usually been considered in media 
with spatially isotropic nonlinear properties. Recently a 
great deal of theoretical and experimental efforts have 
been devoted to studies of nonlinear optical effects and 
soliton formation in photorefractive crystals J^l, 12, 13]. 
While these media exhibit strong nonlinearity at very 
low optical power their nonlinear response is inherently 
anisotropic jl^. The anisotropy causes a number of ob- 
servable effects including astigmatic self-focusing of opti- 
cal beams 15], elliptically shaped solitary solutions 16], 
geometry-sensitive interactions of solitons , and fixed 
optical pattern orientation |l^ . 

Several previous studies of MI in the context of pho- 
torefractive media were limited to a 1-dimensional geom- 
etry where the anisotropy is absent , and the 
physics is similar to the standard saturable nonlinear- 
ity [2^. On the other hand, in a real physical situation 
where one deals with finite sized beams, the anisotropic 
aspects of the photorefractive nonlinear response are ex- 
pected to play a significant role. Some previous work 
[25 [2^ already indicated the importance of anisotropy 
in the transversal break-up of broad beams propagat- 
ing in biased photorefractive crystals. However, no de- 
tailed analysis of this phenomenon was carried out. In 



this paper we study the MI of optical beams in photore- 
fractive media taking into account the full 2-dimensional 
anisotropic model of the photorefractive nonlinearity. 

Time independent propagation of an optical beam 
£{r^z) = (A/2)e*^^^~'^*^ + c.c. in a nonlinear medium 
with a weakly varying index of refraction is governed by 
the parabolic equation 
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Here r = (x^y) and z are transverse and axial coordi- 
nates, = x{d/dx) + y{d/dy)^ k = 27rno/A, A is the 
wavelength in vacuum, uj = 2i\cj\ c is the speed of light, 
and n = no + n2(r, z) is the refractive index, with no the 
spatially uniform background index, and n2 the spatially 
varying nonlinear increment. 

In the case of a photorefractive screening nonlinear- 
ity the optical beam propagates through a photorefrac- 
tive crystal externally biased with a DC electric field. 
The beam excites charges which after migrating due to 
diffusion and drift in the applied field, are subsequently 
trapped by impurity or defect centers. The effective non- 
linearity (refractive index change) is proportional to the 
low frequency electric field E created by light induced 
charge redistribution Q • Iii the situation of in- 

terest here where the optical field is linearly polarized 
along X which coincides with the crystalline c axis the 
nonlinear increment to the refractive index is given by 
^2(^5 z) = — \nQr^^Ex{r^ z), with raa the relevant compo- 
nent of the electro-optic tensor, and Ex the x component 
of the low frequency electric field in the medium. 

It is convenient to describe the nonlinear material re- 
sponse in terms of the quasi-static potential induced by 
the optical field. As shown in the appendix the resulting 
set of dimensionless equations is 
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where / = (l + |Ap)/[l-^V±-(enV±0)]. The coordinates 
and variables have been normahzed using the scahngs 
given in the appendix with the addition of|Ap//^ |Ap 
with / = 2Is/{eQnQc). 

Equation Q describes the most general situation when 
the electrostatic potential in the crystal is induced by 
two distinct transport mechanisms: drift of charges in 
the biasing DC field plus photogalvanic field and their 
diffusion. The relative strength of the diffusion and drift 
terms is determined by the dimensionless parameter 



a = 



The diffusion contribution which leads to spatially asym- 
metric stimulated scattering dominates at large trans- 
verse wavenumbers of order k]j. On the other hand the 
drift terms give the dominant contribution to the spa- 
tially symmetric MI which is prominent at much smaller 
transverse wavenumbers. Thus the term proportional to 
a in Eq. ^ is often neglected when studying MI. 

The initial linear stage of the filamentation instability 
may be investigated by putting 
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where F represents the growth rate of the perturbation 
characterized by the transverse wavevector q, and fre- 
quency detuning f2. 

The steady state solution to Eq. © in the one di- 
mensional plane wave limit with a = = is d<p/dx = 
(1 + Eph/Ee^t)\Ao\'^/{l + l^oP) where Aq is the ampli- 
tude of the plane wave which is assumed to vanish for 
X ±oo. It is therefore convenient to renormalize the 
potential as = (1 + Epi,/Ee^t)\Ao\^x/{l + |^oP) + <^ so 
that the equations of motion become 
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with / = (1 + - • (enV±(^)]. The ansatz (gj 

then gives + ^ph/^ext)|^o| V(l + which 

provides a continuous transition to the ID solution. 

Solving the linear problem defined by Eqs. (|4I6|) gives 
the dispersion relation 
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where en{0) = cos^(^) + sin^(^) with ey the static di- 
electric tensor component along y and the angle of q 
with respect to the x axis. In the limit of a single trans- 
verse dimension {qy = 0) without diffusion Eq.JT)) reduces 
to the formula for the growth rate in saturable nonlinear 
media 

The instability growth rate is given by Re[r] = 
yRelr^JT^/V^. In FiglHwe show the growth rate 
as a function of the spatial frequency q = |q| for a few 
values of the angle using parameters characteristic of 
a photorefractive crystal as given in the appendix. The 
growth rate depends strongly on the angular orientation 
of the initial perturbation. In particular, it always attains 
the largest value when the wave vector of the perturba- 
tion coincides with the direction of the applied electric 
field = 0). As 6> departs from zero the amplification of 
the perturbation decreases, and the growth rate becomes 
less strongly peaked at small q. The growth rate is an 
even function of q provided 1] = and inspection of Eq. 
JZj) shows that it falls of for large q as 1/q. Interestingly 
Eq.Q predicts there is no instability for perturbations 
with wave vectors perpendicular to the direction of the 
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FIG. 1: Growth rate as a function of spatial frequency q for 
few values of the transverse angle for = 1, = 0, 
^ph = 0, ty/tc = 0.53, and ^ = 0.42. The solid lines show 
the diffusionless case {a = 0) and the dashed lines correspond 
to a = 0.046. 



applied field. 

When the perturbation is frequency shifted with re- 
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FIG. 2: Growth rate of frequency shifted perturbations vs 
spatial frequency q for few values of the of transverse angle 
for Aq = 1, Qr = 1, 0, e^/cc 0.53, a 0.046, and 

i = 0.42. 

spect to the plane wave {Vt ^ 0) the growth rate becomes 
an asymmetric function of q. This is shown in Fig. [21 for 
VLt = 1. Positive q in the figure corresponds to the direc- 
tion of a plane wave that experiences two-wave mixing 
gain. We see that the growth rate has several maxima as 
a function of q in the presence of a frequency shift. 

Finally we note that the angular dependence takes on 
a simple form in the small q limit where the diffusion 
contribution can be neglected. Putting a = ^ = in {TJ 
and assuming no frequency shift so ^7 = we obtain 

(8) 

This latter expression shows clearly the decline of the 
instability growth rate with increasing angle away from 
the X axis. 



Ao=0.6 

• 


Ao=0.8 


Ao=1 

'II' 


Ao=1.2 








% 


% 


q. 






w • 









FIG. 3: Intensity distribution (top) and its spatial spectrum 
(bottom)of a gaussian beam with initial random noise, after 
propagation over a distance of 5mm in a photorefractive crys- 
tal. The size of the computational window is 200/i x 200/i. A 
zero-frequency component has been removed from the spec- 
trum. 

In order to verify the reliability of the linear approx- 
imation discussed above we resorted to numerical anal- 
ysis of the full 2-dimensional model governing propaga- 
tion of optical beams in a photorefractive medium Eqs. 



FIG. 4: Intensity distribution (a) and its spatial spectrum (b) 
of a broad Gaussian beam with a superimposed initial ran- 
dom noise, after propagation over a distance of 5mm in an 
isotropic nonlinear medium. For clarity of presentation the 
zero-frequency component has been removed from the spec- 
trum. 



(|2I3|) . For simplicity and comparison with the experi- 
mental conditions discussed below we only considered the 
frequency degenerate {ft = 0) and short Debye length 
{a = ^ = 0) limit. We also assumed the lack of a pho- 
togalvanic effect {Eph = 0). We used a split step Fast 
Fourier transform code to solve the propagation equa- 
tion ((2j and a finite difference technique to find the elec- 
trostatic potential and refractive index change induced 
by the beam. Results of these calculations are shown 
in Fig. (|3l5|) . In all cases the direction of the applied 
DC field is horizontal (along the x-axis). In Fig.Q we 
show a few examples of numerical simulations depicting 
results of propagation of the wide Gaussian beam with 
initial random perturbation of its amplitude, for a few 
values of the beams peak intensity. These graphs illus- 
trate the inherently anisotropic nature of the instability. 
Initially random perturbations lead to amplification of 
the perturbation with almost zero y-component of the 
wave vector. This leads to appearance of beam modu- 
lation in the form of roughly vertically oriented stripes. 
Graphs in the bottom row display the spatial spectrum 
(a zero frequency component has been removed for clar- 
ity of presentation) of the resulting intensity distribution. 
Notice that the spatial frequency with the highest growth 
rate varies with the peak intensity of the beam, which is 
in agreement with the prediction of Eq. (|5|). For compar- 
ison FigEI shows the intensity distribution obtained with 
the Gaussian beam propagating in a standard isotropic 
saturable medium. The beam experiences modulational 
instability but this time all amplified spatial frequencies 
are located on a ring reflecting the isotropic nature of the 
nonlinear process. 

The complete absence of instability for = 7r/2 and 
its apparent one-dimensional character as depicted in 
Fig.© are a direct consequence of the anisotropy of 
the nonlinear response of the photorefractive medium. 
The light induced focusing power is roughly three times 
stronger in the direction of the applied DC field than 
in the direction perpendicular to it p^. Unlike the 
isotropic system where all spatial frequencies correspond- 



ing to highest growth rate are amphfied (Fig^, in the 
photorefractive crystal the highest gain is experienced 
only by perturbations with Qy ^ {6 ^ 0). Therefore 
only these frequencies will contribute to the initial stages 
of the modulational instability described by the linear 
theory. Spatial perturbations with nonzero qy compo- 
nents have much weaker growth rates and will play an 
important role only after the ID structure with Qy ^ 
has reached sufficiently high intensity Then the full 
2-dimensional break-up and subsequent filamentation of 
the beam will follow 23, 25]. However, the full analy- 
sis of such a process is beyond the scope of the present 
paper. 

Next we simulated propagation of a broad Gaussian 
beam with its amplitude perturbed by a spatially pe- 
riodic modulation. The angle which determines the 
angular orientation of the perturbation with respect to 
the direction of applied DC field was varied from ^ = 
to ^ = 90°. The strength of the perturbation (relative 
to the peak intensity of the beam) was less than 10~^. 
Results of the propagation of this beam over a distance 
of five millimetres are shown in Fig|51 Each row of this 
figure corresponds to a different spatial frequency of the 
initial perturbation. A decrease in the amplification of 
the perturbation with increasing angle is evident. 
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FIG. 6: Experimentally observed self-induced modulational 
instability of the initially Gaussian beam (width=200/i). (a) 
light intensity distribution; (b) spatial spectrum (with zero 
frequency component removed). Amplitude of the external 
biasing DC field E^xt — 'ikV/cm. DC field applied in hor- 
izontal direction. The width of the experimental window is 
320/i X 320/i 



applied along the 5 mm long c-axis. The optical beam 
(ImW) from a solid state laser (A = 532 nm) was loosely 
focused at the input face of the crystal. The output inten- 
sity distribution was imaged by a CCD camera and stored 
in a computer. The crystal was illuminated by a broad 
white light beam which was used to control the degree of 
saturation. Typically, the peak intensity of the incident 
beam was of the same order as the average intensity of 
the white light background. We used either an unper- 
turbed beam or beam with superimposed weak periodic 
perturbations. Results of the experiments are shown in 
Fig l6l7l FiglBfa) shows the light intensity distribution at 
the exit facet of the crystal (after 10 mm of propagation) 
in the case where the incoming beam was not intention- 
ally perturbed. It is evident that nonlinear ity induces 
modulational instability which leads to the formation of 
quasi 1-dimensional vertical stripes oriented perpendicu- 
larly to the direction of the applied DC field. Figure|ni(b) 
depicts the corresponding Fourier spectrum of the outgo- 
ing beam where the two distinct peaks with almost zero 
qy components clearly indicate the anisotropic character 
of the instability. The presence of the small y-component 
in the spectrum is the result of a slight misalignment of 
the crystal. 



FIG. 5: Intensity distribution of the optical beam in a biased 
photorefractive crystal with an initially imposed periodic per- 
turbation as a function of the angular orientation of the per- 
turbation and its wave vector {q). Peak intensity DC 
biasing electric field applied along horizontal (x) direction. 
The size of the computational window is 200/i x 200/i. 
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To verify our theoretical findings we conducted exper- 
iments using a crystal of photorefractive strontium bar- 
ium niobate as the nonlinear medium. The experimental 
setup is analogous to that used in our earlier studies of 
photorefractive soliton formation 13]. The crystal was 
5x5x10 mm in size with the optical beam propagating 
along the 10 mm axis and a DC electric field of 1.1 kV 



FIG. 7: Experimentally observed dependence of the modula- 
tional instability of the broad Gaussian beam on the angu- 
lar orientation of initial periodic perturbation. Amplitude of 
the external DC electric field Eext — 2.2kV/cm Field applied 
along the horizontal direction. 

Finally we investigated the role of anisotropy in mod- 



ulational instability of a broad beam with an initially 
imposed periodic perturbation. To this end the incom- 
ing Gaussian beam was initially transmitted through a 
parallel plate which resulted in the appearance of a weak 
spatial periodic modulation of the beam wavefront. By 
rotating the plate we were able to change the orienta- 
tion of this modulation. The perturbed beam subse- 
quently propagated through the biased photorefractive 
crystal. Results of this experiment are shown in Figd 
Grey scale plots in the top row of this figure represent 
the light intensity distribution at the output face of the 
photorefractive crystal corresponding to different angu- 
lar orientations of the periodic pattern characterized by 
the angle {6). Graphs in the bottom row illustrate the 
corresponding intensity profile. As Fig.Q clearly shows 
the amplification of the perturbation decreases rapidly as 
the angle departs from = 0. The rightmost plot shows 
the intensity pattern at the input face of the crystal. For 
better visualisation we plot in Fig. [51 the experimentally 
measured growth rate (normalized to its maximum value) 
as a function of the angle 0. The points represent exper- 
imental data while the line is a theoretical fit (Eq.(|5|) 
with Eph = 0, Ao = 2.6, ^ 0.42, a = and = 0.35. 
Again, the drop in amplification of the perturbation for 
increasing 6 is evident. 
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FIG. 8: Experimentally measured normalized growth rate of 
an initial periodic perturbation imposed onto the broad Gaus- 
sian beam as a function of an angular orientation of the initial 
pattern. Dots - experimental points; line - theoretical fit (for- 
mula (0) 

In conclusion, we investigated modulational instability 
of plane waves and finite beams in photorefractive non- 
linear media biased with a DC electric field. We showed 
that the growth rate of perturbation is affected by the in- 
herent anisotropy of the nonlinear response. It is highest 
for perturbations whose wave- vectors correspond to the 
direction of the biasing DC field. For arbitrarily oriented 
perturbations the effect of anisotropy manifests itself in 
a decrease of the effective strength of the nonlinear re- 
sponse until it reaches zero for wave-vectors perpendic- 
ular to the direction of the field. Our theoretical pre- 



dictions were confirmed by experimental observations in 
strontium barium niobate crystals. 
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APPENDIX A: DERIVATION OF EQUATIONS 

The set of equations describing the optical properties 
of a photorefractive crystal, known as the Kukhtarev 
equations [2^, are 



dt 



^ = {P^CjIem){ND-N^)-7rneN^ (Ala) 



p = e(7V+ -Na- Tie) (Alb) 

(Ale) 

V • (eoeE) = p (Aid) 
^+V-J = 0. (Ale) 

Here Ndj A/"^, Na^ and Ue are the density of donors, ion- 
ized donors, acceptors, and conduction electrons, P and 
a are the coefficients of thermal and photoexcitation, /gm 
is the optical intensity, 7^ is the electron recombination 
coefficient, — e is the charge on an electron, eo is the per- 
meability of vacuum, e is the static dielectric tensor, 
is the Boltzmann constant, T the temperature, p the 
electron mobility, p the charge density, J the current, 
and E the static electric field. Note that the coefficient 
a includes the photogalvanic contribution due to /3ph so 
that we could write the total photoexcitation coefficient 
as cr = ai + Pph. 

We analyze these equations following the approach of 
Ref. 01 . In the absence of thermal or photoexcitation 
ng = so the condition < >=< Na >, where <> 
denotes a spatial average, ensures bulk charge neutrality. 
The negatively charged acceptors do not participate in 
the photoexcitation dynamics so the density Na is fixed 
and serves to limit the magnitude of the photoexcited 
space charge field. To analyze the Kukhtarev equations 
we assume N]j ^ Na ^ ^e- Gauss's law ()Ald|) then 
gives 
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Introducing the Debye wave number for ch arge motion 
along the c axis as /c^ = e^y Na/ (/c^Teoec), where €c is 
the component of the dielectric tensor along c and the 
characteristic field E = ksTko/e = cNa/ {^q^c^d) we 
can write the last expression as 
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with Cn the static dielectric tensor divided by eg. 



The assumption of fast carrier recombination implies 
that due/dt can be set to zero in the equation for charge 
continuity. It fohows that 
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where M = — 1. To proceed we use Eq. ((Ala|) to write 
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where Is = /3/a is the saturation intensity for which the 
rate of thermal excitation equals the rate of photoexcita- 
tion. With fast carrier recombination and Nd ^ Na we 
have 
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with / = lem/h- The characteristic density of electrons 
associated with Ig is no = /^r){ND /Na) 
tron density can be written as 



so the elec- 
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Furthermore the photogalvanic term can be written as 
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Defining the characteristic relaxation time of the elec- 
tric field to 

^p^^-^^ Eq. (El can be written as 



and the photogalvanic field Ep^ 
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vanishes at the boundaries of the medium. Using in- 
stead of E in Eq. Q results in only a small change in 
the wavenumber of the beam which has no physical im- 
portance for this work. We then introduce a potential 
through the relation E^ = —EV^cj) so that Eq. (|A5|) can 
be written as 



to^[V-(enV(/))] + ^ 



M 

^E^dl 

E dxl^ E dx 



1 + M 



^ext S 1+/ 



^V^ft^. (A6) 
kD 1 + M ^ ^ 



In the situation of interest here where the optical 
field is linearly polarized along x the nonlinear incre- 
ment to the refractive index is given by 712 (r, z) = 
— ^nor33£^a;(r, z), with raa the relevant component of the 
electro-optic tensor, and E^ the x component of the 
low frequency electric field in the medium. In a parax- 
ial approximation the optical field therefore satisfies the 
parabolic equation 
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Within the same paraxial approximation we drop the lon- 
gitudinal derivatives in Eq. (|A6|) so that the gradient op- 
erator becomes = xd/dx + yd/dy. Finally introduc- 
ing the spatial scales l± and /y and redefining the coordi- 
nates and variables through {x^y)/l± (^,^), 
(j){E/ EQ^t){l/l±) (/), gives the normalized set 
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This equation coincides with Ref. Eq. (2) with 

X = and 6 = 0. 

We are interested in the situation where the optical 
beam is small compared to the size of the nonlinear 
medium. The externally applied bias field is £^ext = V/ 
with V the applied voltage and the width of the 
medium along x which is taken to coincide with the c 
axis. It is convenient to subtract this field from the op- 
tically induced field so that the field E^ = E — E^^tx 



where r = to//, / 

E,^t/{El±kD), l\\ = 2/(fcngr33^ext), and = ^VM- 
We can calculate characteristic values for the theo- 
retical parameters for the SBN crystal used in the ex- 
perimental work. Using data from Ref. we have 



880, Cy 



470, and Na ^ 10^^ cm^-^^ so that the 



Debye length is 27r/kD = 2.2 jam, and the characteristic 
internal field is ^ = 730 V/cm. The optical parameters 
are no = 2.3, A = .532 /im, while the electro-optic coeffi- 
cient in our crystal was measured to be rss ~ 180 pm/V. 
With a typical applied field of /^ext = 2.2 kV/cm we have 
/|l 350 /im, l± = 2.5 /im, a 0.046, and ^ = 0.42. 
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